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Abstract. We present two types of meta-algorithm that can greatly 
improve the accuracy of existing algorithms for integrating the equations 
of motion of dynamical systems. The first meta-algorithm takes an inte- 
grator that is time-symmetric only for constant time steps, and ensures 
time-symmetry even in the case of varying time steps. The second meta- 
algorithm can be applied to any self-starting integration scheme to create 
time-symmetry, both for constant and for variable time steps, even if the 
original scheme was not time-symmetric. Our meta-algorithms are most 
effective for Hamilton systems or systems with periodic solutions. If the 
system is not Hamiltonian (for example, if some dissipative force exists), 
our methods are still useful so long as the dissipation is small. 



1. Introduction 

Many problems in computational physics are governed by underlying equations 
that are intrinsically time-symmetric. In particular, for any simulation in Hamil- 
tonian dynamics, we can run the movie of our computation equally well forwards 
as backwards, and in both cases obtain a physically allowable solution. Clearly, 
it is desirable to use an integration algorithm that reflects this time-symmetry 
as a built-in property. In that case, intuitively speaking, particles can no longer 
fly 'off the tracks', so to speak, when moving through a curve. An example 
of particular interest to this conference is astrophysical particle simulations, in 
which particles may correspond directly to physical units, such as molecules or 
stars, or may form tracers used to approximate the solution of a set of underlying 
equations with continuous variables. 

The basic idea is this: if an algorithm would make a particle spiral out 
systematically in a given situation, it would have to do so equally in the forward 
and backward direction. Time reversal, however, would force an inward motion, 
and the conclusion is that the net spiral out (or spiral in) has to be zero, at 
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least when averaged over a full period in a periodic system. The error in energy 
thus has to be periodic as well, and cannot build up from orbit to orbit. Of 
course, errors are still made for finite integration step size, but they show up as 
errors in phase. In many applications, phase errors are preferable over errors in 
quantities such as the total energy, that should be conserved. 

Most algorithms that are in common use for orbit integration in simula- 
tions are not time-symmetric. Even those that arc, typically lose their symme- 
try property as soon as one allows variable time steps. In the remainder of this 
paper, we offer two types of meta-algorithm for constructing a larger class of 
time-symmetric algorithms. In §§2,3 wc show how to restore time symmetry for 
algorithms that were symmetric, but lose that property when one uses variable 
time steps. In §§4,5 we show how to create time symmetry even for those algo- 
rithms that were never time symmetric to begin with, not even in the constant 
time step case. 

2. Building A Better Leapfrog 

A celebrated example of an integration scheme with built-in time symmetry 
is the leapfrog scheme, also known as the Verlet method. It is widely used in 
many applications of particle simulations, such as in molecular dynamics, plasma 
physics, fluid dynamics and stellar dynamics (Hockney & Eastwood 1988; Barnes 
&; Hut 1986). The time symmetry is manifest in the interleaved representation, 
which gave rise to the 'leapfrog' name: 



where r can stand for the position vector of a single particle or the combined 
vector ri,r2, . . . , representing a system of N particles. The quantity v = 
dr/dt is the velocity and a{t) = a{r{t)) = dv/dt the acceleration. The subscripts 
after the various quantities indicate the time at which they apply, in units of 
the time step, i. e. vi = v{t + \5t). 

It is convenient to map the standard interleaved description into a form in 
which all variables are defined at the same instant in time: 



Starting from {ro,uo,ao}) one first computes ri, then ai(ri), by evaluating the 
appropriate expression dictated by the system under consideration, and finally 
Vi. While Eq.(2) looks as though it has lost its explicit time symmetry, it is 
still equivalent to the original Eq.(l), as can be verified by direct substitution 
(Barnes & Hut 1989). However, if the time step 6t is allowed to vary, through a 
functional dependence 5t = h{r, v) for example, time symmetry is lost. 

Time symmetry can be restored if we force the time step to be a symmetric 
function of the begin point and end point of each time step, as was first shown 
by Hut et al. (1995). For example, we can use 




(la) 



{Ih) 



ri = rQ + vo6t + ^ao{St)^, 
^'i = ■"0 + l{ao + ai)St. 



(2a) 
(26) 



St = ^[h{ro,vo) + h{ri,vi)] 



(3) 
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With this choice of time step, the combined set of equations (2,3) has become 
impHcit: in order to determine {ri,vi} from {ro,fo}, we need to know the 
time step size St, which in turn is dependent on {ri,vi}. We can solve this 
problem by starting with 5t = h{rQ,VQ) as a first approximation. This will 
give us approximate values for {ri,vi}, from which we can determine a more 
accurate value for 6t using equation (3). If necessary, this iteration process can 
be repeated several times, but in practice one or two iterations are generally 
sufficient to reach time symmetry to machine accuracy. 



3. Restoring Time Symmetry 

The notion of time symmetry restoration for variable time steps, give above for 
the particular case of the leapfrog scheme, carries over to any class of integration 
schemes that is explicitly time symmetric for constant time steps, as shown by 
Hut et al. (1995). Another example, given in the same paper, concerns the 
following natural fourth-order generalization of the leapfrog scheme: 

ri = ro + livi + vo)dt - ^^(ai - ao){6tf, (4a) 

vi = vo + i(«i + ao)St - ^{ji - jo){Stf, (46) 

which is a truncated form of the Hermite scheme (Makino 1991a). Here the jerk 

j = da/dt is calculated directly by differentiation of the expression for the force 
(thereby introducing a dependency on velocity as well as position in the case of 
Newtonian gravitational forces). This set of equations is manifestly time sym- 
metric for constant time steps. Unlike the original second-order leapfrog, given 
in eq. (2), the scheme in eq. (4) is implicit, already for constant time steps. 
Applying the same symmetrization procedure given in equation (3) leads, after 
iteration, to a fully time-symmetric fourth-order generalization of the leapfrog. 
In practice, iteration can give a quick convergence, leading to substantially im- 
proved accuracy for a fixed amount of computer time (Hut et al. 1995). 



4. Creating Time Symmetry 

Even if no time symmetry is present in a given algorithm, it is possible to 
construct an implicit version that is manifestly time symmetric, for the general 
class of self-starting (i. e. one-step) integration schemes. The mcta-algorithm 
that performs this feat is described in detail by Funato et al. (1997), as a 
generalization of the prescription offered by Hut et al. (1995). Here we give a 
brief outline of the main idea behind the treatment. 
Let us start with the ordinary differential equation 

and a given self-starting integration scheme, expressed as 

Vi+i = yi + F{xi,yi;hi), (2) 
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where Xj and yi are the values of the variables x and y after the i-th step, and 
hi is the step size at Xj. 

Our meta-algorithm can be expressed as follows: 



where hi can be constructed in the form of a function hi = /(/ij, /ij+i) that 
is symmetric in its arguments: f{x,y) = f{y,x). For example, we could take 
simply 



an implicit formula for yj+i. As before, we can solve this equation by iteration, 
starting with the original non-symmetric scheme as the initial trial function. 

While this recipe is surprisingly simple, we can do even better. Instead of 
taking a given estimate for the step size, we can use the difference 



to estimate the local truncation error. This information can be used to imple- 
ment a form of adaptive step size control. See Funato et al. (1997) for further 
details. 

5. Building A Better Runge-Kutta 

As an example application of our more general meta-algorithm, we have con- 
structed a time-symmetric version of the popular fourth-order Runge-Kutta in- 
tegration scheme. Figure || shows the behavior of the errors in the relative energy 
and angular momentum for a binary orbit with initial eccentricity e = 0.9 (the 
values plotted are determined at apocenter). The dashed and solid curves show 
the time evolution of the errors for the standard Runge-Kutta scheme and for 
the symmetrized Runge-Kutta scheme, respectively. In both cases, variable step 
sizes have been used. The number of time steps is comparable in both cases 
(around 600 per orbit). 

Figure |^ shows that no discernible secular error is produced, even after 1000 
orbital periods, for the run integrated by the time-symmetric fourth-order Runge 
Kutta Method. In contrast, the error increases linearly for the run integrated 
by the standard fourth-order Runge Kutta Method. Further quantitative details 
will be provided in the cost /performance analysis by Funato et al. (1997). 

6. Discussion 

We have reviewed two types of meta-algorithm, based upon a time-symme- 
trization procedure. The first meta-algorithm preserves time symmetry that 
would otherwise be lost when integration step sizes are allowed to change during 
integration. The second meta-algorithm creates time symmetry, even for those 
algorithms where no symmetry was present in the equal-step-size case. 



Vi+i = 



yi + F{xi,yi]hi) = yi + - F{xi,yi;hi) - F{xi+i,yi+i; -hi) , (3) 





(5) 
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Figure 1. Effects of our meta-algorithm applied to a fourth-order 
Runge-Kutta scheme. Plotted is the growth of the relative error in 
energy (left) and angular momentum (right) for a Kcplcr orbit with 
eccentricity 0.9. Dashed and full lines correspond to the standard 
Runge-Kutta scheme and the time-symmetrized Runge-Kutta version, 
respectively. 



We mention here briefly a few recent applications of these ideas. McMillan 
8z Hut (1996) have constructed a fully automated package for performing grav- 
itational three-body scattering experiments, where the central orbit integrator 
is built along the principles outlined by Hut et al. (1995). They found that 
the symmetrization meta-algorithm gave a significant speed-up to the fourth- 
order Hermite scheme used. Most importantly, they found that the fraction of 
rejected experiments was diminished greatly, compared to the standard integra- 
tion scheme. The problem here is that some resonant scattering experiments 
can stay in an intermediate state for a very long time, before finally decaying 
into a final state. No matter how accurate an initial time step criterion has 
been chosen, there will always been a small fraction of such lingering states that 
will ultimately lead to an unacceptable build-up of errors (both integration er- 
rors and round-off errors). Time symmetrization, while not circumventing this 
problem, can greatly alleviate the situation. 

Another application has been discussed by Funato et al. (1996a,b). They 
have applied time symmetrization to the Kustaanheimo-Stiefel regularization 
method, a sophisticated way to 'unfold' the singularity of the three-dimensional 
Kepler problem by mapping each point in three dimensional space to a unit circle 
in an auxiliary four dimensional space. The combination of these two powerful 
techniques has resulted in the most accurate way yet designed to integrate orbits 
near collision singularities. 

For some applications, specific adaptations of our general meta-algorithm 
can make good use of the known constraints inherent in the underlying problem. 
One example that we have recently explored is that large-scale simulations of 
planetary formation. The problem is that close encounters and physical colli- 
sions of planetesimals make it absolutely necessary to use individually variable 
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time steps (Aarseth 1985). All standard choices for highly accurate integration 
schemes, such as the symplectic schemes, lose their desirable properties once we 
allow individual particles to change their integration time step length at will. In 
contrast, our meta- algorithm shows a way out, as demonstrated by Kokubo &; 
Makino (1997). They made use of the fact that most planctcsimals have nearly 
circular orbits (e << 0.01), which means that their time step is practically con- 
stant when block time step (McMillan 1986, Makino 1991b) are used, even if 
we allow the use of variable time steps. Even though the time step size of a 
particle shrinks significantly during a close encounter, leading to a break-down 
of time-symmetry, this occurs only for a very small fraction of time, for a typical 
particle. Since most of the integration error actually comes from the integration 
of nearly unperturbed orbits around the sun, reserving strict time symmetry for 
unperturbed orbits turns out to be a good compromise, leading to high overall 
accuracy. Details will be provided by Kokubo & Makino (1997). 
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